Capitulo_5
Gabriel Omar Maya Briseño
2025-04-29
5 Practical Issues in the Specification and Estimation of Discrete Choice Models
library(knitr)
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(discrtr) # A companion package for the book Introduction to Discrete Choice Analysis with `R`
library(dplyr) # A Grammar of Data Manipulation
library(ggplot2) # Create Elegant Data Visualisations Using the Grammar of Graphics
library(mlogit) # Multinomial Logit Models## Cargando paquete requerido: dfidx
##
## Adjuntando el paquete: 'dfidx'
##
## The following object is masked from 'package:stats':
##
## filterlibrary(readr) # Read Rectangular Text Data
#library(stargazer) # Well-Formatted Regression and Summary Statistics Tables
library(gplots) # Various R Programming Tools for Plotting Data##
## Adjuntando el paquete: 'gplots'
##
## The following object is masked from 'package:stats':
##
## lowesslibrary(ggmosaic) # Mosaic Plots in the 'ggplot2' Framework
library(treemapify)
library(ggridges)
library(ggalluvial)
library(evd)
library(htmlwidgets) # HTML Widgets for R
library(kableExtra) # Construct Complex Table with kable and Pipe Syntax##
## Adjuntando el paquete: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rowslibrary(plotly) # Create interactive web graphics##
## Adjuntando el paquete: 'plotly'
##
## The following object is masked from 'package:ggplot2':
##
## last_plot
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following object is masked from 'package:graphics':
##
## layoutlibrary(tidyr) # Tidy messy data
#library(webshot2) # Take screenshots of web pagesFunciones de utilidad
data("mc_commute_wide", package = "discrtr")
mc_commute_wide[1:6, 1:10]## id choice available.Cycle available.Walk available.HSR available.Car
## 1 566872636 HSR No Yes Yes No
## 2 566873140 HSR No Yes Yes Yes
## 3 566874266 HSR No No Yes Yes
## 4 566874842 Walk Yes Yes Yes No
## 5 566881170 Walk Yes Yes Yes No
## 6 566907438 Walk No Yes Yes No
## time.Cycle time.Walk time.HSR time.Car
## 1 NA 21.31439 5 NA
## 2 NA 12.78863 10 2
## 3 NA NA 15 4
## 4 5.828157 20.00000 8 NA
## 5 5.828157 20.00000 5 NA
## 6 NA 10.00000 3 NASe toma a continuacion choice como columna y tiempo las variables que se pueden elegir
example_wide <- mc_commute_wide %>% dplyr::select(id, choice, starts_with("time")) |>
filter(id %in% c(566910139, 566873140, 566872636))
example_wide## id choice time.Cycle time.Walk time.HSR time.Car
## 1 566872636 HSR NA 21.31439 5 NA
## 2 566873140 HSR NA 12.78863 10 2
## 3 566910139 Walk 4.371118 15.00000 20 5Usaremos la misma selección de variables y filas para ver la diferencia con una tabla en formato por longitud.
example_wide |> pivot_longer(cols = starts_with("time."), names_prefix = "time.", names_to="alternative", values_to="time")## # A tibble: 12 × 4
## id choice alternative time
## <dbl> <fct> <chr> <dbl>
## 1 566872636 HSR Cycle NA
## 2 566872636 HSR Walk 21.3
## 3 566872636 HSR HSR 5
## 4 566872636 HSR Car NA
## 5 566873140 HSR Cycle NA
## 6 566873140 HSR Walk 12.8
## 7 566873140 HSR HSR 10
## 8 566873140 HSR Car 2
## 9 566910139 Walk Cycle 4.37
## 10 566910139 Walk Walk 15
## 11 566910139 Walk HSR 20
## 12 566910139 Walk Car 5Nótese como el mismo id se repite en cuatro filas, pero el tiempo está dado en una sola. El package mlogit utiliza tablas en formato por longitud. Dado que las tablas en formato por amplitud son más comunes, este package incluye una función para cambiar la tabla. Para ello es necesario indicar cuáles son las variables que varían en cada alternativa, en este caso son cuatro variables: time, access (el tiempo necesario para llegar a una parada de autobús), wait (tiempo de espera por un autobús) y transfer (el número de transferencias cuando se utiliza el autobús). Las últimas tres variables son específicas de HSR (cuando se utiliza el transporte público)
example_long <- mc_commute_wide %>% filter(id %in% c(566910139, 566873140, 566872636)) |>
mlogit.data(shape="wide", choice="choice", varying=3:22)
example_long## ~~~~~~~
## first 10 observations out of 12
## ~~~~~~~
## id choice parking vehind gender age shared
## 1 566872636 FALSE No No Man 21 Living in Shared Accommodations
## 2 566872636 FALSE No No Man 21 Living in Shared Accommodations
## 3 566872636 TRUE No No Man 21 Living in Shared Accommodations
## 4 566872636 FALSE No No Man 21 Living in Shared Accommodations
## 5 566873140 FALSE No Yes Man 23 No
## 6 566873140 FALSE No Yes Man 23 No
## 7 566873140 TRUE No Yes Man 23 No
## 8 566873140 FALSE No Yes Man 23 No
## 9 566910139 FALSE No Yes Man 49 No
## 10 566910139 FALSE No Yes Man 49 No
## family child street_density sidewalk_density LAT LONG
## 1 No No 14.37621 22.63322 43.26302 -79.90074
## 2 No No 14.37621 22.63322 43.26302 -79.90074
## 3 No No 14.37621 22.63322 43.26302 -79.90074
## 4 No No 14.37621 22.63322 43.26302 -79.90074
## 5 No No 19.49754 39.64003 43.25885 -79.90476
## 6 No No 19.49754 39.64003 43.25885 -79.90476
## 7 No No 19.49754 39.64003 43.25885 -79.90476
## 8 No No 19.49754 39.64003 43.25885 -79.90476
## 9 Living with Family Yes 15.39380 26.06793 43.25597 -79.91260
## 10 Living with Family Yes 15.39380 26.06793 43.25597 -79.91260
## PersonalVehComf_SD PersonalVehComf_D PersonalVehComf_A PersonalVehComf_SA
## 1 0 0 1 0
## 2 0 0 1 0
## 3 0 0 1 0
## 4 0 0 1 0
## 5 0 0 0 0
## 6 0 0 0 0
## 7 0 0 0 0
## 8 0 0 0 0
## 9 0 0 1 0
## 10 0 0 1 0
## Fun_SD Fun_D Fun_A Fun_SA ActiveNeigh_SD ActiveNeigh_D ActiveNeigh_A
## 1 0 1 0 0 0 0 1
## 2 0 1 0 0 0 0 1
## 3 0 1 0 0 0 0 1
## 4 0 1 0 0 0 0 1
## 5 0 0 0 1 0 0 0
## 6 0 0 0 1 0 0 0
## 7 0 0 0 1 0 0 0
## 8 0 0 0 1 0 0 0
## 9 0 0 1 0 0 0 0
## 10 0 0 1 0 0 0 0
## ActiveNeigh_SA UsefulTrans_SD UsefulTrans_D UsefulTrans_A UsefulTrans_SA
## 1 0 0 0 1 0
## 2 0 0 0 1 0
## 3 0 0 0 1 0
## 4 0 0 0 1 0
## 5 0 1 0 0 0
## 6 0 1 0 0 0
## 7 0 1 0 0 0
## 8 0 1 0 0 0
## 9 0 0 0 1 0
## 10 0 0 0 1 0
## BusComf_SD BusComf_D BusComf_A BusComf_SA TravelAlone_SD TravelAlone_D
## 1 0 1 0 0 1 0
## 2 0 1 0 0 1 0
## 3 0 1 0 0 1 0
## 4 0 1 0 0 1 0
## 5 0 0 0 0 1 0
## 6 0 0 0 0 1 0
## 7 0 0 0 0 1 0
## 8 0 0 0 0 1 0
## 9 0 0 1 0 0 1
## 10 0 0 1 0 0 1
## TravelAlone_A TravelAlone_SA Shelters_SD Shelters_D Shelters_A Shelters_SA
## 1 0 0 0 0 1 0
## 2 0 0 0 0 1 0
## 3 0 0 0 0 1 0
## 4 0 0 0 0 1 0
## 5 0 0 0 0 1 0
## 6 0 0 0 0 1 0
## 7 0 0 0 0 1 0
## 8 0 0 0 0 1 0
## 9 0 0 0 0 0 0
## 10 0 0 0 0 0 0
## Community_SD Community_D Community_A Community_SA personal_veh_comfortable
## 1 0 0 0 0 A
## 2 0 0 0 0 A
## 3 0 0 0 0 A
## 4 0 0 0 0 A
## 5 0 1 0 0 N
## 6 0 1 0 0 N
## 7 0 1 0 0 N
## 8 0 1 0 0 N
## 9 0 0 1 0 A
## 10 0 0 1 0 A
## getting_there_fun like_active_neighborhood commute_useful_transition
## 1 D A A
## 2 D A A
## 3 D A A
## 4 D A A
## 5 SA N SD
## 6 SA N SD
## 7 SA N SD
## 8 SA N SD
## 9 A N A
## 10 A N A
## buses_comfortable prefer_travel_alone shelter_good_quality sense_community
## 1 D SD A N
## 2 D SD A N
## 3 D SD A N
## 4 D SD A N
## 5 N SD A D
## 6 N SD A D
## 7 N SD A D
## 8 N SD A D
## 9 A D N A
## 10 A D N A
## numna alt available time access wait transfer chid idx
## 1 2 Car No NA 0 0 0 1 1:Car
## 2 2 Cycle No NA 0 0 0 1 1:ycle
## 3 2 HSR Yes 5.000000 3 15 0 1 1:HSR
## 4 2 Walk Yes 21.314387 0 0 0 1 1:Walk
## 5 3 Car Yes 2.000000 0 0 0 2 2:Car
## 6 3 Cycle No NA 0 0 0 2 2:ycle
## 7 3 HSR Yes 10.000000 4 15 0 2 2:HSR
## 8 3 Walk Yes 12.788632 0 0 0 2 2:Walk
## 9 4 Car Yes 5.000000 0 0 0 3 3:Car
## 10 4 Cycle Yes 4.371118 0 0 0 3 3:ycle
##
## ~~~ indexes ~~~~
## chid alt
## 1 1 Car
## 2 1 Cycle
## 3 1 HSR
## 4 1 Walk
## 5 2 Car
## 6 2 Cycle
## 7 2 HSR
## 8 2 Walk
## 9 3 Car
## 10 3 Cycle
## indexes: 1, 2data.frame(example_long) |> dplyr::select(id, choice, alt, starts_with("time"), idx)## id choice alt time idx
## 1 566872636 FALSE Car NA 1:Car
## 2 566872636 FALSE Cycle NA 1:ycle
## 3 566872636 TRUE HSR 5.000000 1:HSR
## 4 566872636 FALSE Walk 21.314387 1:Walk
## 5 566873140 FALSE Car 2.000000 2:Car
## 6 566873140 FALSE Cycle NA 2:ycle
## 7 566873140 TRUE HSR 10.000000 2:HSR
## 8 566873140 FALSE Walk 12.788632 2:Walk
## 9 566910139 FALSE Car 5.000000 3:Car
## 10 566910139 FALSE Cycle 4.371118 3:ycle
## 11 566910139 FALSE HSR 20.000000 3:HSR
## 12 566910139 TRUE Walk 15.000000 3:WalkDado que hay cuatro alternativas, cada fila corresponde a la situación de elección para una alternativa para un individuo (las tablas en formato por longitud por lo general tienen filas por cada decisor). El conjunto de índices contiene dos variables, la primera identifica al individuo y la segunda la alternativa:
data.frame(example_long$idx)## chid alt
## 1 1 Car
## 2 1 Cycle
## 3 1 HSR
## 4 1 Walk
## 5 2 Car
## 6 2 Cycle
## 7 2 HSR
## 8 2 Walk
## 9 3 Car
## 10 3 Cycle
## 11 3 HSR
## 12 3 WalkEn este caso las variables disponibles para el análisis son:
colnames(example_long)## [1] "id" "choice"
## [3] "parking" "vehind"
## [5] "gender" "age"
## [7] "shared" "family"
## [9] "child" "street_density"
## [11] "sidewalk_density" "LAT"
## [13] "LONG" "PersonalVehComf_SD"
## [15] "PersonalVehComf_D" "PersonalVehComf_A"
## [17] "PersonalVehComf_SA" "Fun_SD"
## [19] "Fun_D" "Fun_A"
## [21] "Fun_SA" "ActiveNeigh_SD"
## [23] "ActiveNeigh_D" "ActiveNeigh_A"
## [25] "ActiveNeigh_SA" "UsefulTrans_SD"
## [27] "UsefulTrans_D" "UsefulTrans_A"
## [29] "UsefulTrans_SA" "BusComf_SD"
## [31] "BusComf_D" "BusComf_A"
## [33] "BusComf_SA" "TravelAlone_SD"
## [35] "TravelAlone_D" "TravelAlone_A"
## [37] "TravelAlone_SA" "Shelters_SD"
## [39] "Shelters_D" "Shelters_A"
## [41] "Shelters_SA" "Community_SD"
## [43] "Community_D" "Community_A"
## [45] "Community_SA" "personal_veh_comfortable"
## [47] "getting_there_fun" "like_active_neighborhood"
## [49] "commute_useful_transition" "buses_comfortable"
## [51] "prefer_travel_alone" "shelter_good_quality"
## [53] "sense_community" "numna"
## [55] "alt" "available"
## [57] "time" "access"
## [59] "wait" "transfer"
## [61] "chid" "idx"Iniciaremos definiendo una fórmula que sólo considera el tiempo de transporte.
# Function `mFormula()` is used to define multi-part formulas of the form:
# y ~ x | z | w, which in the notation used for the anatomy of utility functions is
# choice ~ alternative vars. with generic coefficients |
# individual vars. with specific coefficients |
# alternative vars. with specific coefficients
# In this formula time is one of x variables
f1 <- mFormula(choice ~ time)La función model.matrix puede ser usada para examinar cómo la fórmula es aplicada a los datos:
f1 |> model.matrix(example_long)## (Intercept):Cycle (Intercept):HSR (Intercept):Walk time
## 1 0 1 0 5.000000
## 2 0 0 1 21.314387
## 3 0 0 0 2.000000
## 4 0 1 0 10.000000
## 5 0 0 1 12.788632
## 6 0 0 0 5.000000
## 7 1 0 0 4.371118
## 8 0 1 0 20.000000
## 9 0 0 1 15.000000Definamos una fórmula con una variable específica de cada individuo
# Function `mFormula()` is used to define multi-part formulas of the form:
# y ~ x | z | w, which in the notation used for the anatomy of utility functions is
# choice ~ alternative vars. with generic coefficients |
# individual vars. with specific coefficients |
# alternative vars. with specific coefficients
# In this formula `time` is one of x variables and `sidewalk_density` is one of z variables
f2 <- mFormula(choice ~ time | sidewalk_density)La matriz del modelo está dada por el tiempo que se toma usar lod diferentes medios de transporte:
f2 |> model.matrix(example_long)## (Intercept):Cycle (Intercept):HSR (Intercept):Walk time
## 1 0 1 0 5.000000
## 2 0 0 1 21.314387
## 3 0 0 0 2.000000
## 4 0 1 0 10.000000
## 5 0 0 1 12.788632
## 6 0 0 0 5.000000
## 7 1 0 0 4.371118
## 8 0 1 0 20.000000
## 9 0 0 1 15.000000
## sidewalk_density:Cycle sidewalk_density:HSR sidewalk_density:Walk
## 1 0.00000 22.63322 0.00000
## 2 0.00000 0.00000 22.63322
## 3 0.00000 0.00000 0.00000
## 4 0.00000 39.64003 0.00000
## 5 0.00000 0.00000 39.64003
## 6 0.00000 0.00000 0.00000
## 7 26.06793 0.00000 0.00000
## 8 0.00000 26.06793 0.00000
## 9 0.00000 0.00000 26.06793Ahora consideramos que el tiempo es específico de cada alternativa, en lugar de tener coeficientes genéricos, por medio de la fórmula:
f3 <- mFormula(choice ~ 0 | sidewalk_density | time)Nótese que, dado que no definimos variables específicas de cada alternativa con coeficientes genéricos, tenemos que establecerlas como cero. La matriz del modelo es:
f3 %>% model.matrix(example_long)## (Intercept):Cycle (Intercept):HSR (Intercept):Walk sidewalk_density:Cycle
## 1 0 1 0 0.00000
## 2 0 0 1 0.00000
## 3 0 0 0 0.00000
## 4 0 1 0 0.00000
## 5 0 0 1 0.00000
## 6 0 0 0 0.00000
## 7 1 0 0 26.06793
## 8 0 1 0 0.00000
## 9 0 0 1 0.00000
## sidewalk_density:HSR sidewalk_density:Walk time:Car time:Cycle time:HSR
## 1 22.63322 0.00000 0 0.000000 5
## 2 0.00000 22.63322 0 0.000000 0
## 3 0.00000 0.00000 2 0.000000 0
## 4 39.64003 0.00000 0 0.000000 10
## 5 0.00000 39.64003 0 0.000000 0
## 6 0.00000 0.00000 5 0.000000 0
## 7 0.00000 0.00000 0 4.371118 0
## 8 26.06793 0.00000 0 0.000000 20
## 9 0.00000 26.06793 0 0.000000 0
## time:Walk
## 1 0.00000
## 2 21.31439
## 3 0.00000
## 4 0.00000
## 5 12.78863
## 6 0.00000
## 7 0.00000
## 8 0.00000
## 9 15.00000Experimento 1 Consiste en que los valores mu y beta son iguales
ts <- tibble(Individual=1:6, Choice=c("A","A", "B", "A", "B", "B"),
yA=c(1,1,0,1,0,0), yB=c(0, 0, 1, 0, 1, 1),
xA=c(5,2,5,1,4,3), xB=c(4, 5, 2, 6, 1, 4))
# Set the parameters:
mu <- 0
beta <- 0
my_prob <- function(xA, xB){
exp(beta * xA)/(exp(beta * xA) + exp(mu + beta * xB))
}
# Calculate probabilities. Notice that these are the logit probabilities
# Individual 1
P1A <- my_prob(ts$xA[1], ts$xB[1])
P1B <- 1 - P1A
# Individual 2
P2A <- my_prob(ts$xA[2], ts$xB[2])
P2B <- 1 - P2A
# Individual 3
P3A <- my_prob(ts$xA[3], ts$xB[3])
P3B <- 1 - P3A
# Individual 4
P4A <- my_prob(ts$xA[4], ts$xB[4])
P4B <- 1 - P4A
# Individual 5
P5A <- my_prob(ts$xA[5], ts$xB[5])
P5B <- 1 - P5A
# Individual 6
P6A <- my_prob(ts$xA[6], ts$xB[6])
P6B <- 1 - P6A
# Calculate likelihood function as the product of all the probabilities
# Each probability is raised to ynj
L <- P1A^ts$yA[1] * P1B^ts$yB[1] *
P2A^ts$yA[2] * P2B^ts$yB[2] *
P3A^ts$yA[3] * P3B^ts$yB[3] *
P4A^ts$yA[4] * P4B^ts$yB[4] *
P5A^ts$yA[5] * P5B^ts$yB[5] *
P6A^ts$yA[6] * P6B^ts$yB[6]
# Create data frame to tabulate results:
df_experiment_1 <- data.frame(Individual = c(1, 2, 3, 4, 5, 6),
Choice = c("A", "A", "B", "A", "B", "B"),
PA = c(P1A, P2A, P3A, P4A, P5A, P6A),
PB = c(P1B, P2B, P3B, P4B, P5B, P6B))
# Display table
kable(df_experiment_1, "html", digits = 4, booktabs = TRUE, align = c("c", "c", "c", "c")) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
footnote(general = paste("The value of the likelihood function in Example 1 is: ",
round(L, digits = 4)))| Individual | Choice | PA | PB |
|---|---|---|---|
| 1 | A | 0.5 | 0.5 |
| 2 | A | 0.5 | 0.5 |
| 3 | B | 0.5 | 0.5 |
| 4 | A | 0.5 | 0.5 |
| 5 | B | 0.5 | 0.5 |
| 6 | B | 0.5 | 0.5 |
| Note: | |||
| The value of the likelihood function in Example 1 is: 0.0156 |
Experimento 2 1 Consiste en que los valores mu y beta son diferentes
# Set the parameters:
mu <- 0.5
beta <- -0.5
my_prob <- function(xA, xB){
exp(beta * xA)/(exp(beta * xA) + exp(mu + beta * xB))
}
# Calculate probabilities. Notice that these are the logit probabilities
# Individual 1
P1A <- my_prob(ts$xA[1], ts$xB[1])
P1B <- 1 - P1A
# Individual 2
P2A <- my_prob(ts$xA[2], ts$xB[2])
P2B <- 1 - P2A
# Individual 3
P3A <- my_prob(ts$xA[3], ts$xB[3])
P3B <- 1 - P3A
# Individual 4
P4A <- my_prob(ts$xA[4], ts$xB[4])
P4B <- 1 - P4A
# Individual 5
P5A <- my_prob(ts$xA[5], ts$xB[5])
P5B <- 1 - P5A
# Individual 6
P6A <- my_prob(ts$xA[6], ts$xB[6])
P6B <- 1 - P6A
# Calculate likelihood function as the product of all the probabilities
# Each probability is raised to ynj
L <- P1A^ts$yA[1] * P1B^ts$yB[1] *
P2A^ts$yA[2] * P2B^ts$yB[2] *
P3A^ts$yA[3] * P3B^ts$yB[3] *
P4A^ts$yA[4] * P4B^ts$yB[4] *
P5A^ts$yA[5] * P5B^ts$yB[5] *
P6A^ts$yA[6] * P6B^ts$yB[6]
# Create data frame to tabulate results:
df_experiment_2 <- data.frame(Individual = c(1, 2, 3, 4, 5, 6),
Choice = c("A", "A", "B", "A", "B", "B"),
PA = c(P1A, P2A, P3A, P4A, P5A, P6A),
PB = c(P1B, P2B, P3B, P4B, P5B, P6B))
# Display table
kable(df_experiment_2, "html", digits = 4, booktabs = TRUE, align = c("c", "c", "c", "c")) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
footnote(general = paste("The value of the likelihood function in Example 2 is: ",
round(L, digits = 4)))| Individual | Choice | PA | PB |
|---|---|---|---|
| 1 | A | 0.2689 | 0.7311 |
| 2 | A | 0.7311 | 0.2689 |
| 3 | B | 0.1192 | 0.8808 |
| 4 | A | 0.8808 | 0.1192 |
| 5 | B | 0.1192 | 0.8808 |
| 6 | B | 0.5000 | 0.5000 |
| Note: | |||
| The value of the likelihood function in Example 2 is: 0.0672 |
Al cambiar los valores provoca que cambien las probabilidades
# Create a grid to plot the likelihood function
mu = seq(from = -1, to = 1, by = 0.01)
beta = seq(from = -2, to = 0, by = 0.01)
coeffs <- expand.grid(mu, beta)
my_prob <- function(xA, xB, mu, beta){
exp(beta * xA)/(exp(beta * xA) + exp(mu + beta * xB))
}
ts <- data.frame(Individual = c(1, 2, 3, 4, 5, 6),
Choice = c("A", "A", "B", "A", "B", "B"),
yA = c(1, 1, 0, 1, 0, 0),
yB = c(0, 0, 1, 0, 1, 1),
xA = c(5, 2, 5, 1, 4, 3),
xB = c(4, 5, 2, 6, 1, 4))
# Define the likelihood function
lkh <- function(mu = 0, beta = 0){
P1A <- my_prob(ts$xA[1], ts$xB[1], mu, beta)
P1B <- 1 - P1A
P2A <- my_prob(ts$xA[2], ts$xB[2], mu, beta)
P2B <- 1 - P2A
P3A <- my_prob(ts$xA[3], ts$xB[3], mu, beta)
P3B <- 1 - P3A
P4A <- my_prob(ts$xA[4], ts$xB[4], mu, beta)
P4B <- 1 - P4A
P5A <- my_prob(ts$xA[5], ts$xB[5], mu, beta)
P5B <- 1 - P5A
P6A <- my_prob(ts$xA[6], ts$xB[6], mu, beta)
P6B <- 1 - P6A
P1A^ts$yA[1] * P1B^ts$yB[1] *
P2A^ts$yA[2] * P2B^ts$yB[2] *
P3A^ts$yA[3] * P3B^ts$yB[3] *
P4A^ts$yA[4] * P4B^ts$yB[4] *
P5A^ts$yA[5] * P5B^ts$yB[5] *
P6A^ts$yA[6] * P6B^ts$yB[6]
}
# Evaluate the likelihood function on the grid
L <- lkh(mu = coeffs$Var1, beta = coeffs$Var2)
L <- data.frame(mu = coeffs$Var1, beta = coeffs$Var2, L)
L <- xtabs(L ~ beta + mu, L) %>% # Convert to cross-tabulation matrix
unclass() # Drop the xtabs class (plotly does not like it)
likelihood_plot <- plot_ly(z = ~L, x = ~mu, y = ~beta) %>%
add_surface() %>%
layout(scene = list(
xaxis = list(title = "x-axis (mu)"),
yaxis = list(title = "y-axis (beta)"),
zaxis = list(title = "z-axis (L)")))
likelihood_plotEn la gráfica se puede observar que los valores aproximados que maximizan la función de verosimilitud son mu = 0.10 y bet = 0.65. Si utilizamos estos coeficientes para calcular las probabilidades, podemos comparar las probabilidades de los experimentos 1 y 2:
# Set the parameters:
mu <- 0.1
beta <- -0.65
my_prob <- function(xA, xB){
exp(beta * xA)/(exp(beta * xA) + exp(mu + beta * xB))
}
# Calculate probabilities. Notice that these are the logit probabilities
# Individual 1
P1A <- my_prob(ts$xA[1], ts$xB[1])
P1B <- 1 - P1A
# Individual 2
P2A <- my_prob(ts$xA[2], ts$xB[2])
P2B <- 1 - P2A
# Individual 3
P3A <- my_prob(ts$xA[3], ts$xB[3])
P3B <- 1 - P3A
# Individual 4
P4A <- my_prob(ts$xA[4], ts$xB[4])
P4B <- 1 - P4A
# Individual 5
P5A <- my_prob(ts$xA[5], ts$xB[5])
P5B <- 1 - P5A
# Individual 6
P6A <- my_prob(ts$xA[6], ts$xB[6])
P6B <- 1 - P6A
# Calculate likelihood function as the product of all the probabilities
# Each probability is raised to ynj
L <- P1A^ts$yA[1] * P1B^ts$yB[1] *
P2A^ts$yA[2] * P2B^ts$yB[2] *
P3A^ts$yA[3] * P3B^ts$yB[3] *
P4A^ts$yA[4] * P4B^ts$yB[4] *
P5A^ts$yA[5] * P5B^ts$yB[5] *
P6A^ts$yA[6] * P6B^ts$yB[6]
# Create data frame to tabulate results:
df_approx_solution <- data.frame(Individual = c(1, 2, 3, 4, 5, 6),
Choice = c("A", "A", "B", "A", "B", "B"),
PA = c(P1A, P2A, P3A, P4A, P5A, P6A),
PB = c(P1B, P2B, P3B, P4B, P5B, P6B))
# Join tables for displaying results
df <- df_experiment_1 %>% left_join(df_experiment_2,
by = c("Individual", "Choice")) %>%
left_join(df_approx_solution,
by = c("Individual", "Choice"))
# Display table
kable(df,"html", digits = 4, booktabs = TRUE, col.names = c("Individual", "Choice", "PA", "PB", "PA", "PB", "PA", "PB"),
align = c("l", "c", "c", "c", "c", "c", "c", "c")) %>%
kable_styling(latex_options = c("striped")) %>%
add_header_above(c(" " = 1, " " = 1,
"Experiment 1" = 2,
"Experiment 2" = 2,
"Approx Max Likelihood" = 2))%>%
footnote(general = paste("The approximate optimal value of the likelihood function is: ",
round(L, digits = 4)))| Individual | Choice | PA | PB | PA | PB | PA | PB |
|---|---|---|---|---|---|---|---|
| 1 | A | 0.5 | 0.5 | 0.2689 | 0.7311 | 0.3208 | 0.6792 |
| 2 | A | 0.5 | 0.5 | 0.7311 | 0.2689 | 0.8641 | 0.1359 |
| 3 | B | 0.5 | 0.5 | 0.1192 | 0.8808 | 0.1141 | 0.8859 |
| 4 | A | 0.5 | 0.5 | 0.8808 | 0.1192 | 0.9589 | 0.0411 |
| 5 | B | 0.5 | 0.5 | 0.1192 | 0.8808 | 0.1141 | 0.8859 |
| 6 | B | 0.5 | 0.5 | 0.5000 | 0.5000 | 0.6341 | 0.3659 |
| Note: | |||||||
| The approximate optimal value of the likelihood function is: 0.0763 |
Dado que la función de verosimilitud toma valores entre 0 y 1,la función log-likelihood toma valores menores que cero.
modelo logit
Consideraremos nuevamente el modelo para elección de modo de transporte. Es necesario transformar todo el conjunto de datos del formato por amplitud por longitud:
mc_commute_long <- mc_commute_wide |>
mlogit.data(shape="wide",
# Name of column with the choices
choice = "choice",
# Numbers of columns with attributes that vary by alternative
varying = 3:22)Esta función requiere al menos dos argumentos: una mFormula y un conjunto de datos. Podemos verificar que las fórmulas definidas previamente son objetos de este tipo:
class(f1)## [1] "mFormula" "Formula" "formula"class(f2)## [1] "mFormula" "Formula" "formula"class(f3)## [1] "mFormula" "Formula" "formula"# Function `mlogit()` is used to estimate logit models
# It needs a multi-part formula and a data set in long form
model1 <- mlogit(f1, mc_commute_long)
# Function `summary()` give the summary of data objects,
# including the output of model estimation algorithms
summary(model1)##
## Call:
## mlogit(formula = choice ~ time, data = mc_commute_long, method = "nr")
##
## Frequencies of alternatives:choice
## Car Cycle HSR Walk
## 0.204364 0.034182 0.244364 0.517091
##
## nr method
## 5 iterations, 0h:0m:0s
## g'(-H)^-1g = 0.000102
## successive function values within tolerance limits
##
## Coefficients :
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept):Cycle 0.366898 0.201434 1.8214 0.06854 .
## (Intercept):HSR 0.707267 0.118104 5.9885 2.118e-09 ***
## (Intercept):Walk 3.212834 0.184169 17.4451 < 2.2e-16 ***
## time -0.056494 0.005685 -9.9374 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -768.63
## McFadden R^2: 0.50323
## Likelihood ratio test : chisq = 1557.2 (p.value = < 2.22e-16)Los resultados incluyen la proporción de frecuencia observada de cada alternativa además de información acerca del procedimiento de optimización.
Ahora se estima otro modelo usando la segunda fórmula:
model2 <- mlogit(f2, mc_commute_long)
summary(model2)##
## Call:
## mlogit(formula = choice ~ time | sidewalk_density, data = mc_commute_long,
## method = "nr")
##
## Frequencies of alternatives:choice
## Car Cycle HSR Walk
## 0.204364 0.034182 0.244364 0.517091
##
## nr method
## 5 iterations, 0h:0m:0s
## g'(-H)^-1g = 0.000287
## successive function values within tolerance limits
##
## Coefficients :
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept):Cycle -0.3058620 0.4541379 -0.6735 0.500629
## (Intercept):HSR 0.6761001 0.2150745 3.1436 0.001669 **
## (Intercept):Walk 2.3453233 0.3124706 7.5057 6.106e-14 ***
## time -0.0562191 0.0056910 -9.8785 < 2.2e-16 ***
## sidewalk_density:Cycle 0.0273085 0.0169331 1.6127 0.106803
## sidewalk_density:HSR 0.0016950 0.0085746 0.1977 0.843299
## sidewalk_density:Walk 0.0350392 0.0111749 3.1355 0.001715 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -760.42
## McFadden R^2: 0.50854
## Likelihood ratio test : chisq = 1573.7 (p.value = < 2.22e-16)Es posible seleccionar el nivel de referencia para las utilidades cuando se estima el modelo, en el siguiente ejemplo se utiliza la alternativa “walk” como referencia:
model2 <- mlogit(f2, mc_commute_long, reflevel = "Walk")
summary(model2)##
## Call:
## mlogit(formula = choice ~ time | sidewalk_density, data = mc_commute_long,
## reflevel = "Walk", method = "nr")
##
## Frequencies of alternatives:choice
## Walk Car Cycle HSR
## 0.517091 0.204364 0.034182 0.244364
##
## nr method
## 5 iterations, 0h:0m:0s
## g'(-H)^-1g = 0.000287
## successive function values within tolerance limits
##
## Coefficients :
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept):Car -2.3453233 0.3124706 -7.5057 6.106e-14 ***
## (Intercept):Cycle -2.6511853 0.4187557 -6.3311 2.434e-10 ***
## (Intercept):HSR -1.6692232 0.2302925 -7.2483 4.221e-13 ***
## time -0.0562191 0.0056910 -9.8785 < 2.2e-16 ***
## sidewalk_density:Car -0.0350392 0.0111749 -3.1355 0.001715 **
## sidewalk_density:Cycle -0.0077307 0.0148205 -0.5216 0.601935
## sidewalk_density:HSR -0.0333442 0.0084297 -3.9556 7.635e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -760.42
## McFadden R^2: 0.50854
## Likelihood ratio test : chisq = 1573.7 (p.value = < 2.22e-16)Observemos que ahora dos coeficientes de sidewalk son significativos. Mientras la densidad de acera no cambia significativamente la utilidad entre utilizar bicicleta o caminar, vivir en un lugar con alta densidad de aceras reduce la utilidad de usar automóvil y transporte público respecto a caminar.
summary(mc_commute_long$sidewalk_density)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 18.19 22.63 24.18 35.70 59.41Posteriormente, del conjunto de datos utilizados para estimar el modelo seleccionamos algunas observaciones para explorar las densidades de aceras en el rango entre 0 y 60, en intervalos de longitud 5. Por lo tanto, se seleccionan 52 observaciones (filas) de la tabla (trece niveles de densidad de aceras por cuatro alternativas):
mc_commute_predict <- mc_commute_long[1:52, ]Reemplazamos la variable sidewalk density usando valores entre 0 y 60, en intervalos de longitud 5. Dado que cada alternativa es una fila, es necesario crear una sucesión de valores repetidos de la siguiente manera:
mc_commute_predict$sidewalk_density <- rep(seq(from=0, to=60, by=5), each=4)Podemos examinar los valores asignados para sidewalk density
mc_commute_predict |> data.frame() |> select(sidewalk_density) |> slice_head(n=10)## sidewalk_density
## 1 0
## 2 0
## 3 0
## 4 0
## 5 5
## 6 5
## 7 5
## 8 5
## 9 10
## 10 10El conjunto de datos para predicción ahora incluye el rango de valores de sidewalk density en los cuales se tiene interés. El objetivo de la simulación es calcular la probabilidad de escoger alguna de las opciones de transporte conforme la densidad de acera varía, para ello también es necesario fijar los valores de otras variables, en particular se fija el valor del tiempo (se considera la mediana):
median(mc_commute_predict$time, na.rm=TRUE)## [1] 10mc_commute_predict$time <- 10Ahora se calculan las probabilidades por medio de la función predict y el modelo 2:
probs <- predict(model2, newdata=mc_commute_predict)El resultado (value) de predict es una matriz que contiene las probabilidades para trece niveles de densidad de acera y cuatro alternativas de transporte:
print(probs)## Walk Car Cycle HSR
## 1 0.7381288 0.07072471 0.05208792 0.13905856
## 2 0.7646590 0.06149222 0.05191414 0.12193469
## 3 0.7887761 0.05323772 0.05152105 0.10646509
## 4 0.8105461 0.04591523 0.05093563 0.09260301
## 5 0.8300785 0.03946495 0.05018525 0.08027133
## 6 0.8475135 0.03381834 0.04929656 0.06937162
## 7 0.8630106 0.02890248 0.04829466 0.05979230
## 8 0.8767385 0.02464350 0.04720262 0.05141540
## 9 0.8888677 0.02096925 0.04604115 0.04412190
## 10 0.8995646 0.01781113 0.04482852 0.03779578
## 11 0.9089873 0.01510533 0.04358057 0.03232680
## 12 0.9172834 0.01279350 0.04231084 0.02761229
## 13 0.9245881 0.01082299 0.04103075 0.02355815Para facilitar el procedimiento al realizar la gráfica, se agregan los valores de la densidad de acera y se cambia la forma de la tabla (por longitud), de tal manera que cada fila sea la probabilidad de la combinación algún modo de transporte y densidad de acera:
probs <- data.frame(sidewalk_density=seq(from=0, to=60, by=5), probs) %>%
pivot_longer(cols=-sidewalk_density,
names_to="Mode",
values_to = "Probability")
probs %>% slice_head(n=10)## # A tibble: 10 × 3
## sidewalk_density Mode Probability
## <dbl> <chr> <dbl>
## 1 0 Walk 0.738
## 2 0 Car 0.0707
## 3 0 Cycle 0.0521
## 4 0 HSR 0.139
## 5 5 Walk 0.765
## 6 5 Car 0.0615
## 7 5 Cycle 0.0519
## 8 5 HSR 0.122
## 9 10 Walk 0.789
## 10 10 Car 0.0532Ahora podemos generar la gráfica
ggplot(probs)+
geom_line(aes(x=sidewalk_density, y=Probability, color=Mode), linewidth=1.2)+
labs(x= expression("Sidewalk density (km/km"^2*")"), y="Probability")
es denominado el modelo “Market Shares”:
f0 <- mFormula(choice ~ 1)
model0 <- mlogit(f0, mc_commute_long)
summary(model0)##
## Call:
## mlogit(formula = choice ~ 1, data = mc_commute_long, method = "nr")
##
## Frequencies of alternatives:choice
## Car Cycle HSR Walk
## 0.204364 0.034182 0.244364 0.517091
##
## nr method
## 6 iterations, 0h:0m:0s
## g'(-H)^-1g = 9.89E-06
## successive function values within tolerance limits
##
## Coefficients :
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept):Cycle -1.788207 0.157592 -11.3470 < 2e-16 ***
## (Intercept):HSR 0.178756 0.080839 2.2113 0.02702 *
## (Intercept):Walk 0.928318 0.070464 13.1743 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -1547.2
## McFadden R^2: 0
## Likelihood ratio test : chisq = 0 (p.value = 1)La verosimilitud logarítmica del modelo Market Shares es -1547.2, mientras que, para el modelo 2 es -760.418. Luego, la de McFadden para el modelo 2 es:
1-as.numeric(model2$logLik)/as.numeric(model0$logLik)## [1] 0.5085353En el resumen de los modelos, la prueba de proporción de verosimilitud está reportada, y se contrasta con el modelo nulo. Cuando es necesario comparar dos modelos no nulos se utiliza la función lrtest:
lrtest(model1, model2)## Likelihood ratio test
##
## Model 1: choice ~ time
## Model 2: choice ~ time | sidewalk_density
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 4 -768.63
## 2 7 -760.42 3 16.429 0.0009258 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1